Método de la transformada inversa

Daremos por hecho que podemos generar variables aleatorias independientes e identicamente distribuidas con distribución Uniforme$(0,1)$.

Empezamos por la distribución más sencilla posible: $\text{Bernoulli}(p=0.5)$.

Teniendo acceso a una variable aleatoria Uniforme$(0,1)$; ¿cómo simularían la variable aleatoria $X\sim \text{Bernoulli}(0.5)$?

En efecto $X=\mathbb{1}_{\left\{ U \leq 0.5 \right\}}$ se distribuye $\text{Bernoulli}(p=0.5)$

En general si tenemos una distribución discreta $X$ tal que para $x_i\in \left\{ x_1,x_2,\ldots ,x_n\right\}$

\begin{align*} \mathbb{P}\left[X=x_i \right]=p_i \end{align*}

con

\begin{align*} \sum_{i=1}^n p_i =1. \end{align*}

¿Cómo se puede simular $X$ haciendo uso de una variable aleatoria $U\sim\text{Uniforme}(0,1)$?

En efecto podemos definir

\begin{align*} X=\begin{cases} x_1 & \text{ si } U\leq p_1 \\ x_2 & \text{ si } p_1 < U \leq p_1 + p_2 \\ \vdots \\ x_n & \text{ si } p_1 + p_2 + \ldots + p_{n-1}< U \end{cases} \end{align*}

La metodología anterior puede ser generalizada haciendo uso de $F^{\leftarrow}$, la inversa generalizada de la función cumulativa de distribución $F$ asociada a la variable alaeatoria $X$ que se quiere simular.

Recordamos que la función de distribución cumulativa de una variable aleatoria $X$ está dada por

\begin{align*} F(x) = \mathbb{P}\left[ X \leq x \right]. \end{align*}

Una función de distribución cumulativa debe ser no-decreciente, continua por la derecha y tal que

\begin{align*} \lim_{x\to-\infty}F(x)=0 \end{align*}

y \begin{align*} \lim_{x\to \infty}F(x)=1. \end{align*}

Definimos la inversa generalizada de una función cumulativa de distribución como

\begin{align*} F^{\leftarrow}(u)=\inf\left\{ x \, : \, F(x) \geq u \right\} =\min\left\{ x \, : \, F(x) \geq u \right\}, \quad 0<u<1. \end{align*}

Proposición

Dada $F$ una función de distribución cumulativa:

  • a) $u\leq F(x) \iff F^{\leftarrow}(u) \leq x$.
  • b) Si $U\sim\text{Uniforme}(0,1)\implies F^{\leftarrow}(U)$ tiene función cumulativa de distribución $F$ ( es decir $\mathcal{L}(X)=\mathcal{L}(F^{\leftarrow}(U))$).
  • c) Si $F$ es continua $\implies F(X)\sim \text{Uniforme}(0,1)$.

Ejemplo

Si $X\sim \text{Exponencial}(\lambda) \implies F(x)=1-e^{-\lambda x}\implies F^{\leftarrow}(x)=F^{-1}(x)=-\frac{\log(1-x)}{\lambda}$, en pracica se usa $X=-\frac{\log(U)}{\lambda}$ para simular variables aleatorias exponenciales.

In [21]:
using Plots, RandomNumbers.Xorshifts
r = Xoroshiro128Plus(0x1234567890abcdef)
Out[21]:
Xoroshiro128Plus(0xe7eb72d97b4beac6, 0x9b86d56534ba1f9e)
In [22]:
n = 1000
x = -log.(rand(r,n));
In [24]:
plot(p)
Out[24]:
In [26]:
plot(p)
Out[26]:

Simulación de variables aleatorias condicionadas a un intervalo

Si se condiciona $X$ a $X\in (a,b)$ para $a<b$, entonces para $x\leq b$

\begin{align*} \mathbb{P}[X < x | X\in(a,b)] = \frac{ \mathbb{P}[a < X < \min\left\{ x , b\right\} ] }{\mathbb{P}[a<X<b]} = \frac{ F(x) -F(a)}{F(b)-F(a)}. \end{align*}

Se sigue que $F^{-1}(F(a)+U(F(b)-F(a)))\sim X|X\in (a,b)$.

Problema con el método de la transfromación inversa

  • Si $F^{\leftarrow}$ no es analíticamente conocida se tiene que usar aproximaciones a esta función que suelen ser costosas de evaluar. Aún más, la inversa generalizada de $F$ puede tener dependencia compleja de los parámetros en el sentido de que para cierto valor de estos la inversa es menos costosa que para otros.

Método de aceptación rechazo

Sea $X$ una variable aleatoria continua, con función de densidad de probabilidad $f(x)$, de la cual deseamos generar simulaciones. La idea del método de aceptación-rechaza es usar una variable aleatoria auxiliar $Y$, de la cual podemos simular y un criterio de aceptación de la variable aleatoria tal que la distribución de $Y$ condicionada al criterio de aceptación coincida con la distribución de interés asociada a $X$.

En particular consideraremos $Y$ una variable aleatoria auxiliar continua con densidad de probabilidad $g$ y una constante $0<M\in \mathbb{R}$ tales que $f(x) \leq Mg(x) \, \forall x\in \mathbb{R}$.

Algoritmo de Aceptación-Rechazo

1) Simule $Y\sim \mathcal{L}(g)$ y $U\sim \text{Uniforme}(0,1)$.

2) Acepte $X=Y$ si se cumple que $$U\leq \frac{f(Y)}{Mg(Y)};$$ en caso contrario regrese a 1).

Proposición

  • a) Sea $A$ el evento de aceptar en algoritmo de aceptación-rechazo, entonces $$\mathbb{P}\left[ A \right]=\frac{1}{M}.$$
  • b) Si $X$ está generado por el algoritmo de aceptación-rechazo, entonces $$\mathbb{P}[X\leq x]=\int_{-\infty}^x f(z)\mathrm{d}z=F(x).$$

Intuición usando el teorem fundamental de la simulación

Observemos que una función de densidad de probabilidad $f$ siempre se puede escribir como

\begin{align*} f(x) = \int_0^{f(x)}\mathrm{d}u \end{align*}

por lo que $f$ es la densidad marginal de la distribución conjunta

\begin{align*} (X.U) \sim \text{Uniforme}\left\{ (x,u) \, : \, 0<u<f(x) \right\} \end{align*}

Teorema fundamental de la simulación

$$ X\sim \mathcal{L}(f(x)) $$

es equivalente a $$ (X,U)\sim \text{Uniforme}\left\{ (x,u) \, : \, 0<u<f(x) \right\} $$

Simular una dupla $(X,U)$ como en el teorema fundamental es sencillo si simulamos $X\sim \mathcal{L}(f)$ y $U | X \sim \text{Uniforme}(0,f(X))$, sin embargo esto no es de utilidad cuando la simulación de interes corresponde a $\mathcal{L})(f)$. Alternativamente, una metodología de aceptación-rechazo interesante en este contexto consiste en considerar $(Y,V)\sim \text{Uniforme}(M)$ con $M$ un conjunto tal que $\left\{ (x,u) \, : \, 0<u<f(x) \right\}\subseteq M$ y aceptar si la dupla $(Y,V)\in \left\{ (x,u) \, : \, 0<u<f(x) \right\}$,

Para ilustrar tal método de aceptación-rechazo consideramos una densidad $f$ tal que $$ \int_a^b f(x)\mathrm{d}x=1 $$ para $a<b$ números reales. Entonces podemos considerar $m=\max\left\{ f(x) \, : \, a \leq x \leq b \right\} $ y $M=\left\{ (y,v) \, : \, 0<v<m \wedge a<y<b \right\}$. Usamos el siguiente algoritmo de simulación.

  • 1) Simular $Y\sim \text{Uniforme}(a,b)$ y $V\sim \text{Uniforme}(0,m) $ para generar una dupla $(Y,V)$
  • 2) Si $V\leq f(Y)$ aceptamos $X=Y$ en caso contrario regresar al paso 1).
In [34]:
plot(p)
Out[34]:
In [37]:
plot(p)
Out[37]:

Anteriormente consideramos el conjunto $M=\left\{ (y,v) \, : \, 0<v<m \wedge a<y<b \right\}$ que nos ayudo a definir variables auxiliares distribuidas uniformemente sobre $M$. Más en general podemos considerar.

$M=\left\{ (y,v) \, : \, 0<v<m(y) \right\}$. Donde sin pérdida de generalidad podemos escribir $m(y)=Mh(y)$ para $0<M<\infty$ y $h$ una densidad de probabilidad.

Observemos que con $V | Y \sim \text{Uniforme}(0,m(Y))$ la regla de aceptación "$ X=Y \iff V\leq f(Y)$" tiene sentido para un algoritmo de simulación de $f$ si $f(x)\leq m(x)\, \forall x\in \mathbb{R}$; aún más en tal caso se tiene que $\mathcal{L}(X)=\mathcal{L}(f) \iff h=g$.

De la observación anterior para $g$ densidad de probabilidad y $M$ constante positiva tales que $f(x)\leq Mg(x) \, \forall x\in \mathbb{R}$ obtendríamos el siguiente algoritmo

1) Simule $Y\sim \mathcal{L}(g)$ y $V\sim \text{Uniforme}(0,Mg(y))$.

2) Acepte $X=Y$ si se cumple que $$V\leq f(Y)$$ en caso contrario regrese a 1).

Que coincide con el algoritmo de aceptación-rechazo al considerar que $U\sim \text{Uniforme}(0,1)\iff KU\sim \text{Uniforme}(0,K).$

In [40]:
plot(p)
Out[40]:

Observación

Para el algoritmo de aceptación-rechazo es necesario poder evaluar $Mg(x)$ y deseable que esta evaluación tenga un bajo costo computacional. Más aún el algoritmo puede ser modificado para basarse en la evaluación de $cf(x)$ para alguna constante $c>0$, esto permite hacer simulación cuando la constante de normalización de la densidad $f$ es inaccessible o costosa de evaluar.

Algoritmo de Aceptación-Rechazo 2

1) Simule $Y\sim \mathcal{L}(g)$ y $U\sim \text{Uniforme}(0,1)$.

2) Acepte $X=Y$ si se cumple que $$U\leq \frac{cf(Y)}{Mg(Y)};$$ en caso contrario regrese a 1).

En este caso la probabilidad de aceptar es $\frac{c}{M}$ y se mantiene que $X\sim \mathcal{L}(f)$.

Podemos considerar

\begin{align*} f(x) & \propto e^{-x^2/2}\left( sin(6x)^2 +3cos(x)^2 sin(4x)^2 + 1 \right) \\ g(x) & \propto e^{-x^2/2} \end{align*}

y $M=5$

In [42]:
plot(p)
Out[42]:
In [43]:
p=histogram(y_vec[accept_ind], nbins=250, normalize=true, label="")
title!("Histograma de simulaciones de f con aceptación rechazo")
Out[43]:
In [44]:
plot(p)
Out[44]:

Usar una función para "ensandwichar" en aceptación-rechazo

Si podemos evaluar de manera no costosa una función $g_l$ tal que

$$ g_l(x) \leq f(x) \leq M g_m(x) $$

con $M>0$ constante positiva y $g_m$ densidad de probabilidad tales que se puede implementar el algoritmo de aceptación rechazo. Entonce tiene sentido proponer el siguiente algoritmo cuando evaluar $g_l$ es menos costoso computacionalmente que evaluar $f(x)$.

Algoritmo de Aceptación-Rechazo-Sandwich

1) Simule $Y\sim \mathcal{L}(g_m)$ y $U\sim \text{Uniforme}(0,1)$.

2) Acepte $X=Y$ si se cumple que $$U\leq \frac{g_l(Y)}{Mg_m(Y)};$$ en caso contrario proceda a 3).

3) Acepte $X=Y$ si se cumple que $$U\leq \frac{f(Y)}{Mg_m(Y)};$$ en caso contrario regrese a 1).

Densidades Log-cóncavas

Decimos que una función de de densidad de probabilidad es log-cóncava si $h(x)=log(f(x))$ es cóncava, si $f$ es suficientemente suave esto es equivalente a que $h'(x)=\frac{\mathrm{d}h}{\mathrm{d}x}$ sea estrictamente decreciente o $h''(x)=\frac{\mathrm{d}^2 h}{\mathrm{d}x^2}<0$. Para este tipo de densidades es sencillo construir funciones que "ensandwichen" a la densidad de interés.

In [247]:
plot(p)
Out[247]:
In [249]:
plot(p)
Out[249]:
In [276]:
plot(p)
Out[276]:

Recordamos que si $h$ es diferenciable la recta tangente a la gráfica de $h$ en $x_i$ está dado por $T(x)=h(x_i)+h'(x_i)(x-x_i)$.

Dados $x_i<x_{i+1}$, las rectas tangentes a $h$ en $x_i$ y $x_{i+1}$ se intersectan en el punto dado por

\begin{align*} z_i = \frac{h(x_{i+1}) -h(x_i) -x_{i+1}h'(x_{i+1}) +x_ih'(x_i) }{h'(x_i)-h'(x_{i+1})} \end{align*}

Por lo que dados puntos $T_n = \left\{ x_1, \ldots , x_n \right\}$ definimos para $x\in(z_{i-1},z_i]$

$$ M g_u(x) = h(x_i) + h'(x_i)(x-x_i) $$

con $i\in \left\{1, \ldots, n\right\}$, $z_0$ el mínimo del soporte de $h$, en caso de que $h$ tenga soporte acotado por debajo ó $-\infty$ en caso contrario, y $z_k$ es el máximo del soporte de $h$ en caso de que $h$ tenga soporte acotado por arriba ó $\infty$ en caso contrario.

Por otro lado definimos para $x\in [x_i,x_{i+1}]$

$$ g_l(x) = \frac{(x_{j+1}-x)h(x_j)+(x-x_j)h(x_{j+1})}{x_{j+1}-x_j} $$

y si $x<x_1$ o $x>x_n$ definimos $g_l(x)=-\infty$.

Algoritmo de Aceptación Rechazo Adaptativo

1) Incialize $n$ y $T_n$ (el conjunto puntos nudo).

2) Genere $X\sim \mathcal{L}\left( g_u(x) \right)$ y $U \sim \text{Uniforme}(0,1)$.

3) Si $$U\leq \frac{ g_l(X) }{Mg_u(X)}$$ acepte $X$; en caso contrario si $$U\leq \frac{ f(X) }{Mg_u(X)} $$ acepte $X$ y actualize $n=n+1$, $T_{n+1}=T_n \cup \left\{ X \right\}$; en caso contrario regrese a 2).

In [255]:
@time prueb_prev = adaptive_accept_reject_vprev([-1.1,1.1],10000,h,dh);
 29.384700 seconds (1.61 G allocations: 26.664 GiB, 14.47% gc time)
In [256]:
@time prueb = adaptive_accept_reject([-1.1,1.1],10000,h,dh);
  0.606691 seconds (12.94 M allocations: 319.221 MiB, 12.00% gc time)
In [269]:
@time randn(10000);
  0.000094 seconds (6 allocations: 78.359 KiB)
In [292]:
plot(p)
Out[292]:
In [295]:
plot(p)
Out[295]:
In [299]:
plot(p)
Out[299]:

Ejemplo

La distribución gaussiana inversa generalizada está determinada por la siguiente función de densidad de probabilidad

\begin{align*} f_{GIG}(x)= \frac{\left(a/b\right)^{p/2}}{2\mathcal{K}_p\left(\sqrt{ab}\right)} x^{p-1}e^{-\left(ax +b/x \right)/2}\mathbb{1}_{\left\{ x >0 \right\}} \end{align*}

con $a,b>0$, $p\in \mathbb{R}$ y $\mathcal{K}_p$ la función modificada de Bessel del segundo tipo.

Observemos que

\begin{align*} f_{GIG}(x) \propto x^{p-1}e^{-\left(ax +b/x \right)/2}\mathbb{1}_{\left\{ x >0 \right\}} \end{align*}

por lo que existe $c>0$ tal que \begin{align*} h(x)=\log \left( c f_{GIG}(x) \right) =(p-1)\log(x)-\frac{ax}{2} - \frac{b}{2x} \end{align*} satisface \begin{align*} h''(x)=-\frac{(p-1)}{x^2} -\frac{b}{x^3} \end{align*} por lo que $h$ es cóncava para $p\geq1$ con $x$ restringida a $(0,\infty)$.

In [301]:
plot(p)
Out[301]:
In [308]:
plot(p)
Out[308]:

Ejemplo

Regresión Poisson

Supongamos que tenemos muestras $(Y_1,X_1), \ldots , (Y_n,X_n)$ tales que

\begin{align*} Y_i | a,b,X_i & \sim \text{Poisson}\left( e^{aX_i + b} \right) \\ a & \sim \text{Normal}(0,\sigma^2) \\ b & \sim \text{Normal}(0,\tau^2) \end{align*}

Si consideramos observaciones $Y_i=y_i$, la distribución posterior de $a,b|\pmb{Y}=\pmb{y},\pmb{X}=\pmb{x}$ entonces tiene densidad dada por

\begin{align*} \pi(a,b|\pmb{y},\pmb{x}) &\propto \pi(\pmb{y}|\pmb{x},a,b) \pi(a) \pi(b) \\ & \propto exp\left( a\sum_{i=1}^n y_i + b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} -\frac{b^2}{2\tau^2} \right) \end{align*}

Entonces

\begin{align*} log\left( \pi(a|\pmb{y},\pmb{x},b) \right) & = a\sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} \\ log\left( \pi(b|\pmb{y},\pmb{x},a) \right) & = b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{b^2}{2\tau^2} \end{align*}

Podemos considerar \begin{align*} h_1(a) = & a\sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} \\ h_2(b)= & b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{b^2}{2\tau^2} \end{align*} con lo que \begin{align*} h_1'(a) & = \sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a}{\sigma^2} \\ h_2'(b) & = \sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n x_i e^{x_i b} -\frac{b}{\tau^2} \end{align*}

In [6]:
using RDatasets
prussian = dataset("pscl", "prussian")
┌ Info: Recompiling stale cache file /home/workstation/.julia/compiled/v1.2/RDatasets/JyIbx.ji for RDatasets [ce6b1742-4840-55fa-b093-852dadbb1d8b]
└ @ Base loading.jl:1240
Out[6]:

280 rows × 3 columns

YYearCorp
Int32Int32Cat…
1075G
2276G
3277G
4178G
5079G
6080G
7181G
8182G
9083G
10384G
11085G
12286G
13187G
14088G
15089G
16190G
17091G
18192G
19093G
20194G
21075I
22076I
23077I
24278I
25079I
26380I
27081I
28282I
29083I
30084I
In [189]:
# Año
x = sort(unique(prussian[!,:Year]))
Out[189]:
20-element Array{Int32,1}:
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
In [190]:
# Número de muertes por patada de caballo
y = [ sum(prussian[!,:Y][ prussian[!,:Year] .== x[i] ]) for i in 1:length(x) ]
Out[190]:
20-element Array{Int64,1}:
  3
  5
  7
  9
 10
 18
  6
 14
 11
  9
  5
 11
 15
  6
 11
 17
 12
 15
  8
  4
In [310]:
plot(p)
Out[310]:
In [306]:
plot(p)
Out[306]:

Simulación de la cola de una distribución

Sea $s\in\mathbb{R}$ y $X\sim \mathcal{L}\left( f(x) \right)$ con $f$ una función de densidad de probabilidad en $\mathbb{R}$ que es monótona en $[s,\infty)$.

Consideremos $U$ una variable aleatoria que toma valores en $(0,1]$, $\sigma >0$ y $Y=g(U)=s-\sigma\log(U)$ entonces $g^{-1}(y)=e^{ (s-y)/\sigma}$ y $U$ tiene densidad

\begin{align*} f_U(u) & = \frac{f_Y(y)}{\left| \frac{\mathrm{d}u}{\mathrm{d}y} \right|} = \frac{f_Y\left( y\right)}{\left| \frac{\mathrm{d}g^{-1}(y)}{\mathrm{d}y} \right|} = \frac{f_Y(y)}{\left| \frac{1}{g'\left( g^{-1}(y) \right) } \right|} = f_Y(y)\left| g'\left( g^{-1}(y) \right) \right| \\ & = f_Y(y) \left| -\frac{\sigma}{g^{-1}(y)} \right| =\sigma f_Y(y)e^{ (y-s)/ \sigma} = \frac{ \sigma f_Y(y) }{ u } \end{align*}

Por construcción la variable aleatoria $Y$ toma valores en $[s,\infty )$. Por lo que podemos determinar la ley de probabilidad de $U$ al elegir $f_Y(y)=f$, nuestra densidad de probabilidad de interés. Así las cosas, si $(Z_1,Z_2)$ es tomado uniformemente debajo de la gráfica de $f_U(u)$ entonces $Z_1$ sería una realización de $\mathcal{L}(f)$ restringida a $[s,\infty )$, es decir usamos aceptación-rechazo. Este método es especialmente conveniente si con $f_Y=f$ se tiene que $f_U$ es estrictamente decreciente en $[s,\infty)$ ya que entonces podemos usar la cota $Mg(x)=f_U(s)$.

Ejemplo

Si $f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$ entonces $f_U(u)=\frac{ \sigma }{\sqrt{2\pi}}e^{-y^2/2}e^{(y-s)/\sigma} $ es estrictamente decreciente si y sólo si

\begin{align*} \frac{\mathrm{d}}{\mathrm{d}y}\left( e^{-y^2/2 + (y-s)/\sigma} \right) <0 & \iff e^{-y^2/2 + (y-s)/\sigma}\left( -y + 1/\sigma \right) < 0 \\ & \iff \frac{1}{y} < \sigma. \end{align*}

Recordemos que $s < y$ por lo que se puede elegir $\sigma\in [1/s, \infty )$. Así las cosas podemos utilizar $Mg(u)=\sigma e^{-s^2/2}/\sqrt{2\pi}$ y la regla de aceptación

\begin{align*} U_2 \leq \frac{f_U(U_1)}{f_U(s)}=\frac{ e^{-Y^2/2 + s^2/2} }{U_1}=e^{-Y^2/2 + s^2/2 + (Y-s)/\sigma}. \end{align*}

La probabilidad de aceptación en este caso es un factor de $1/\sigma$ al ser $\sigma$ factor de $Mg$. Por esta razón tomamos el mínimo valor posible de $\sigma$ que es $1/s$.

Simulación de una Normal restringida a la cola $[s,\infty )$

1) Genera $U_1, U_2 \sim \text{Uniforme}(0,1)$ independientes y $Y=s - \frac{1}{s}\log(U_1)$.

2) Si $$U_2 \leq e^{-Y^2/2 + s^2/2 +s(Y-s) }$$ aceptamos $X=Y$, en caso contrario regresamos a 1).

Observando que $$U_2 \leq e^{-Y^2/2 + s^2/2 +s(Y-s) } \iff -2 log(U_2) \leq \left( \frac{\log(U_1)}{s} \right)^2$$ obtenemos un algoritmo menos costoso:

1) Genera $U_1, U_2 \sim \text{Uniforme}(0,1)$ independientes y $Y= -\frac{1}{s}\log(U_1)$.

2) Si $$-2\log(U_2) > Y^2$$ aceptamos $X=s+Y$, en caso contrario regresamos a 1).